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Abstract: A finite element method (FEM) for solving the complex valued 
k(a)) vs. (0 dispersion curve of a 3D metamaterial/photonic crystal system 
is presented. This 3D method is a generalization of a previously reported 2D 
eigenvalue method |[Tl|2l. This method is particularly convenient for analyz- 
ing periodic systems containing dispersive (e.g., plasmonic) materials, for 
computing isofrequency surfaces in the k-space, and for calculating the de- 
cay length of the evanescent waves. Two specific examples are considered: 
a photonic crystal comprised of dielectric spheres and a plasmonic fishnet 
structure. Hybridization and avoided crossings between Mie resonances 
and propagating modes are numerically demonstrated. Negative index 
propagation of four electromagnetic modes distinguished by their symmetry 
is predicted for the plasmonic fishnets. By calculating the isofrequency con- 
tours, we also demonstrate that the fishnet structure is a hyperbolic medium. 
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1. Introduction 

The numerical simulation of electromagnetic fields inside both metamaterial and photonic crys- 
tals is an important tool for analyzing these periodic structures. In particular, the eigenmodes of 
crystals, defined as freely propagating waves not coupled to external currents, are often of the 
most interest. The conventional method 13] ID of numerically solving for crystal eigenmodes 
is to define the geometry of the unit cell of the crystal of interest and the differential equation 
that the fields must obey in this geometry and then impose Bloch periodic boundary conditions. 
The partial differential equation problem is then discretized using one of the many standard 
methods (finite element, finite integral, finite difference, etc.) thereby turning it into an alge- 
braic eigenvalue problem with a finite number of degrees of freedom and the frequency co as 
the eigenvalue. This finite sized eigenvalue problem is then solved numerically. An important 
detail of this method is that the Bloch wavenumber k is chosen beforehand, and the frequency is 
then computed as a function of the wavenumber, yielding the dispersion curves o = a)(k). This 
is the most commonly used method for calculating dispersion curves of the electromagnetic 
waves propagating in photonic crystals or in closely related metamaterial crystals. 

There are however, many instances where it is more convenient to specify the frequency (O 
and solve for the wavenumber as a function of frequency: k = k(a)). At least four such instances 
can be identified. First, metamaterials often contain dispersive materials such as metals, where 
the dielectric function strongly depends on the frequency on the wave. In this case, the eigenfre- 
quency problem becomes a nonlinear eigenvalue problem and must be solved iteratively 15]. In 



contrast, when solving for the wavenumber as a function of frequency, the resulting eigenvalue 
problem only needs to be solved once. Second, it is often useful to solve for the wavenumber 
as the eigenvalue because of the information contained in the complex wavenumber including 
the decay lengths of the electromagnetic modes (either due to finite dissipation or because of 
the evanescent nature of the mode) and the figure of merit ||6] of negative index modes. Third, 
in the majority of experiments electromagnetic fields inside metamaterial/photonic crystals are 
excited by external sources producing time-harmonic fields with real valued frequencies. A 
complex wavenumber eigenvalue simulation provides the correct field distribution in the pho- 
tonic crystal relevant to such an experiment. Fourth, this approach provides a very natural way 
of calculating the so-called isofrequency surfaces corresponding to (B(k) = const, where (O is 
real. Isofrequency diagrams are fundamentally important for predicting wave refraction at the 
PC interfaces Q and for calculating the density of states (jS). 

There are several previously published methods on calculating k(G)) dispersion curves in- 
cluding variations of the plane-wave expansion method ||9] [Tol and diagonalizing the crystal 
transfer matrix ifTTI [T2l . A method of solving for complex wavenumber dispersion curves us- 
ing the FEM has been proposed fT] |2l but only for 2D crystals. The benefits of the complex 
wavenumber 2D FEM are becoming better appreciated and use of this method is becoming 
more common lITSl [141 [TSl [T6l [TtI . It is therefore timely to generalize this method of Complex 
Wavenumber Eigenvalue Simulations (CWES) from two to three dimension, which is the object 
of this paper This 3D complex wavenumber eigenvalue simulation was recently used as part of 
a metamaterial homogenization procedure fTSl. 

The basic theory behind solving for complex wavenumber eigenvalues using the FEM dis- 
cretization is explained in Sec. |2] The underlying field equations for the electric and magnetic 
fields and the necessary boundary conditions are discussed. In Sec.[3]we apply this method to a 
3D photonic crystal consisting of non-dispersive dielectric spheres. The photonic band structure 
for electromagnetic waves propagating both parallel and obliquely to the principal axes is cal- 
culated, and different types of modes (transverse and longitudinal) are identified. In Sec. |4]we 
calculate the dispersion curves for a negative index fishnet metamaterial ||T9l and demonstrate 
the existence of four distinct negative index modes. We also calculate the two-dimensional 
isofrequency contours for the least damped transverse mode and demonstrate from the shape of 
the isofrequency contours that this transverse mode is hyperbolic. 

2. The finite element eigenvalue problem 

2.1. The field equation 

In this section we present the FEM formulation for solving for the magnetic field. Electromag- 
netic wave propagation is described by the Maxwell equations which can be rearranged into a 
wave equation for either the electric field E or the magnetic field H. The wave equation for the 
magnetic field is 



Here £(x) and /i(x) are the microscopic permittivity and permeability of the metamate- 
rial/photonic crystal of interest. Due to the periodic nature of the crystal both are are assumed 
to be scalar functions periodic in the crystal lattice. According to Bloch's theorem ||3] SJ the 
magnetic field can be represented as the product of a periodic function and an exponential 
factor 




(1) 



H(x) = u(x) exp[i(c(Jf - k • x)] 



(2) 



where co is the frequency of the wave and k is the wavevector of the Bloch-' 
is a vector function which is periodic in the crystal lattice. By inserting Eq. 
obtain an equivalent field equation for u 



■Floquet wave. u(x) 
^ into Eq. ([T]i we 



^u-^(k-u)-/kx Qvxu^ -/Vx ^\xu^ +Vx Qv x -n^u = 0, (3) 

which can be interpreted as an eigenvalue problem and solved for the Bloch wavenumber k 
as the eigenvalue. The spatial profile of the eigenmode u(x) is also recovered providing the 
magnetic field profile according to Eq. (|2]i and the electric field profile according to 

E(x) = — ^— V X H = — ^— (-ik X u + V X u) exp\i( (Ot - k • x)l (4) 
ieco/c ieco/c 

2.2. The finite element model 

There are several commercial FEM software programs (COMSOL Multiphysics by COMSOL, 
HESS by Ansys, Vector Fields Opera by Cobham Technical Services, etc.) that are available for 
modeling metamaterial/photonic crystals. These commercial software packages provide a con- 
venient graphical user interface for defining a crystal's geometry, meshing the computational 
domain, and visualizing the electromagnetic fields. This allows for models to be quickly devel- 
oped and tested. Of the many commercial FEM codes available, the authors of this paper are 
only aware of one (COMSOL Multiphysics) that allows the user to specify the field equation 
to be solved. The simulation examples and results presented here were obtained using COM- 
SOL. In what follows, only the most essential features of the FEM approach are reviewed; more 
detailed treatments can be found elsewhere ll20l 1211 l22l . 

The FEM is based on setting the integral of a so-called weak expression over the domain 
of interest to zero. Doing so ensures the field equation is satisfied and also creates boundary 
conditions. The weak expression corresponding to Eq. (|3} is 
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Fh(v,u) = -^v -u - - (k- v) (k-u) -i-v • [k X (V X u)] -i (V X v) • - (k X u) 



(5) 
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where v(x) is a test function. When the integral of the weak expression over the unit cell H of 
the crystal is set to zero, integrating by parts gives us two separate integrals (Eq. (|6)). The first 
integral enforces the field equation. The second integral is over the boundary of the domain and 
represents a natural boundary condition | 
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V X -V X u -jU^u 
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n X - (— ik X u + V X u) 
e 



where n is the vector normal to the boundary. On an external boundary, the natural boundary 
condition enforced by the integral in Eq. (|6) over the boundary dD. forces the expression n x 
(— ik X u + V X u) /e to be equal to zero. Recalling Eq. we note that this simply enforces the 
boundary condition n x E = 0. This is known as the perfect electric conductor or PEC boundary 
condition. This natural boundary condition is the default if no other boundary condition is 
enforced. On an internal boundary within the unit cell the surface integrals over each side of 
the boundary must be equal to each other The effect is that the tangential components of the 
electric field must be continuous across the internal boundary or n x = n x where E^ 
and E^ are the electric fields on opposite sides of the internal boundary. 

The periodicity of u is enforced by imposing periodic boundary conditions on the exterior 
boundaries of the unit cell. In COMSOL, these periodic boundary conditions override the nat- 
ural boundary condition. However, if a PEC boundary condition is desired inside the unit cell 
(e.g., on the surface of a metal inclusion) this can be accomplished by removing the subdomain 
representing the metal inclusion. Now only the exterior side of the metal boundary remains and 
the natural boundary condition forces the tangential electric fields to zero on this boundary. 

If a perfect magnetic conductor or PMC boundary condition (n x H) is desired while solving 
for the magnetic field, this can be enforced with constraints |l20l on the tangential magnetic 
field on the boundary. 

In order to turn Eq.(|6} into an eigenvalue problem, the three degrees of freedom that com- 
prise the Bloch wavevector k must be reduced to one by restricting two degrees of freedom. 
This is accomplished by setting k = ko + Ak„ where A will be the eigenvalue solved for, ko 
is an offset vector and k„ is a unit vector (k„ • k„ = 1) that defines the direction of the com- 
plex wavenumber eigenvalue A . The FEM turns the weak form and accompanying boundary 
conditions into an algebraic problem, in this case a quadratic eigenvalue problem 1231 : 

Au + XBu + X^Cu = Q, (7) 

where A, B and C are N x N matrices and m is an N x 1 solution vector N is the number of 
degrees of freedom of the discretized system. Terms in the weak form (Eq. (|5]l) that are zero, 
first and second order in A contribute to the A, B and C matrices respectively. This quadratic 
eigenvalue problem can be recast in the form of a generalized eigenvalue equation 

When using COMSOL to solve the FEM problem, this linearization is performed automatically 
during the solution stage. 

2.3. Electric field formulation 

The previous discussion focused on solving for the magnetic field H or rather the periodic 
function u equal to the magnetic field with the exponential Bloch factor removed. This is espe- 
cially convenient when an inclusion requires a PEC boundary condition since that is the natural 
boundary condition when solving for H. However, solving for the electric field is very similar 
to solving for the magnetic field. The wave equation for the electric field for a free wave is 

Vx QvxE^ -£^E = 0. (9) 
As before, we replace the electric field with a periodic vector field times an exponential factor 



E(x) =u(x)exp[i(a)f-k-x)], 



(10) 



producing the new field equation 



— u--(k-u)-/kx f-Vxu^ -iVx f-kxu^ + Vx f-V x -e^u = 0. (11) 
Mi" VM y VM / VM / c2 

The corresponding weak form for this field equation is 

k'^ I I 1 

Fe(v,u) = — v-u (k • v) (k-u) - i— V- [k X (V X u)] - i(V X v) ■ — (k X u) 



1 aP- 
(V X v) • — (V X u) - e^v-u, 



(12) 



which is equivalent to Eq. (|5]l if e and pL are interchanged. Integrating this weak form over the 
crystal unit cell by parts and setting its value to zero again gives produces two integrals, one 
enforcing the field equation and a surface integral enforcing the boundary condition n x H = 0. 
Thus the PMC boundary condition is the natural boundary condition when solving for the 
electric field. 

3. Example: Dielectric Photonic Crystal 

For a demonstration of the CWES method of calculating complex k dispersion curves we use 
a simple photonic crystal as an example. The unit cell, pictured in Fig. [T] is a cube with a 
dielectric sphere at the center surrounded by vacuum. The sphere has a radius of 0.3fl, where a 
is the lattice constant of the cubic array, and a permittivity of e = 5 — iO.Ol. 

As mentioned in Sec. 12.21 it is necessary to restrict two of the three degrees of freedom of the 
Bloch wavevector k. There are many possible ways to do this. As the first example, we calculate 
the dispersion curves corresponding to the propagation along a principal axis. To simulate this 
we set ko = and k„ = x. The results of this eigenvalue simulation for the frequency range 
\c/a < (0 < 5.5c/a are plotted in Fig.[T]as (0\s. kx=x-ii = X. 

For clarity, we have only plotted the three least evanescent (i.e. possessing the smallest values 
of Im(A:^)) eigenmodes. The three eigenmodes in Fig. [T] are described as either transverse or 
longitudinal according to their polarization. The symmetry of the dispersion curves is such that 
for every solution k(a)) there is also the solution — k(o) indicating that this is a reciprocal 
crystal. The dispersion curves for the transverse modes plotted in Fig. [T] in fact represent two 
polarization-degenerate modes because of the symmetry of the crystal. The longitudinal mode 
with the passband near o = 4.5c/fl is magnetically polarized in the x direction making it a 
magnetic bulk plasmon. The longitudinal mode with the passband near o) = Sc/a is electri- 
cally polarized in the x direction and is an electric bulk plasmon. The field profiles of both 
longitudinal modes indicate that the passbands correspond to Mie's resonances of the dielectric 
sphere ll24l . 

The transverse mode dispersion curve has a band in the approximate frequency range 
4.6c/fl < o < 4-.Ec/a with a large value of lm{kx), indicating it is an evanescent band, but 
a Re(fcT) that is equal to neither nor n/a as is typical of (B(k) dispersion curves. As described 
in Refs. IlllJsl for a quadratic eigenvalue problem with hermitian matrices (corresponding to 
a lossless crystal) the eigenvalues must always be real or come in complex-conjugate pairs. 
The dielectric photonic crystal under consideration has very low loss, so this lossless condi- 
tion approximately holds true for the dispersion curves in Fig. [T] The transverse band in the 
4.6c/fl < © < 4.8c/a frequency band is one half of a complex conjugate pair, the other half is 
a transverse doubly degenerate mode not shown here. At the frequency of o w 4.8c/fl the two 
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Fig. 1. Complex k dispersion curves and field profiles for eigenmodes of the photonic 
crystal pictured in the inset assuming ko = and k„ = x. (a) Real part of A:f(a)) for a 
transversely polarized mode and a diagram of the crystal unit cell, (b) Imaginary part of 
kx{co) for a transversely polarized mode and a field profile for the z polarized transverse 
mode. There are two transverse modes, y and z electrically polarized, which are degenerate, 
(c) Real part of kx{(o) for two longitudinally polarized modes and a field profile for the 
magnetic longitudinal mode, (d) Imaginary part of kx(co) for two longitudinally polarized 
modes and a field profile for the electric longitudinal mode. The longitudinal mode with the 
passband near ft) = A.5c/a is magnetically polarized in the x direction and the longitudinal 
mode with the passband near m = 5c/a is electrically polarized in the x direction. The 
longitudinal modes correspond to Mie's dipole resonances. For all dispersion curves the 
dotted lines are the result of a conventional &)(k) eigenvalue simulation. For all field profiles 
the frequency is CO = 2c /a with arrows representing D,, and D. and color representing D^. 
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Fig. 2. Complex wavenumber dispersion curves and field profiles for eigenmodes of the 
photonic crystal pictured in Fig. [Tla) assuming ko = tt)/csin(;r/6)y and k„ = x. Modes 
excited by p or s polarized incident light are plotted with solid or dashed lines respectively, 
(a) Real part of kx{(o) for two transverse hybrid modes and a expanded view of the avoided 
crossing in Re(fcj;) space, (b) Imaginary part of kx{a)) for two transverse hybrid modes, 
an expanded view of the avoided crossing in \m{kx) space (plotting the same modes as 
the expanded view in Ke.(kx) space), and a field profile for the E- polarized transverse hy- 
brid mode, (c) Real part of kx{c£)) for two longitudinal hybrid modes and a field profile for 
the magnetic longitudinal hybrid mode, (d) Imaginary part of kx{(o) for two longitudinal 
hybrid modes and a field profile for the electric longitudinal hybrid mode. The magnetic 
longitudinal hybrid mode is excited by s polarized incident light and the electric longitudi- 
nal hybrid mode is excited by p polarized incident light. For all field profiles the frequency 
is O) = 2c /a with arrows representing D,. and D- and color representing D^ . 



modes that make up this complex conjugate pair both enter a passband and split, the plotted 
mode going to the F point and the unplotted mode going to the band edge (this unplotted mode 
corresponds to the dotted lines from the ft)(k) simulation). Note that in this passband there are 
two pairs of propagating doubly polarization degenerate modes or four propagating modes in 
total. 

The transverse eigenmodes plotted in Fig.[T]can be excited by a plane wave normally incident 
onto the vacuum/photonic crystal interface if the interface is parallel to the y-z plane. The 
longitudinally polarized modes could not be excited by a normally incident wave without the 
aid of a coupling device at the interface. If the incident beam of light is not normal to the 
interface, if for example the incident beam has a wavenumber laying in the x-y plane but at 
a 30° angle from normal then a different set of eigenmodes will be excited at the interface. 
To simulate these excited eigenmodes we set ko = Co/csm{n/6)y and k„ = x and solve the 
resulting eigenvalue problem. The resulting photonic band structure is plotted in Fig. |2] 

The eigenmodes in Fig. |2] are roughly split into predominantly-transverse (transverse hy- 
brid) and predominantly-longitudinal (longitudinal hybrid) modes. The hybridization between 
the transverse and longitudinal modes is caused by the finite symmetry-breaking Both the 
transverse and longitudinal hybrid modes in Fig. |2]can be characterized by the polarization of 
the incident light that couples to them. For a p polarized incident beam (electric field in the 
x-y plane) the p polarized eigenmode is excited (plotted in Fig. |2] with solid lines) and for an s 
polarized incident beam (electric field in the z direction) the s polarized eigenmode is excited 
(plotted in Fig. |2] with dashed lines). 

At the frequency O) « 4c/a the transverse hybrid modes and the longitudinal hybrids modes 
appear to cross in a propagating band. An expanded view of this region in Fig.|2]plotting both 
transverse and longitudinal hybrid modes shows that the apparent crossing actually occurs in 
a band gap. Viewed in complex space it is clear that this is actually an avoided crossing 
and in the band gap the transverse and longitudinal hybrid modes form complex conjugate 
pairs mull]. 

We see that even for a simple photonic crystal the CWES method of calculating the dispersion 
curve produces rich and complex results. In particular, it is not possible to solve for evanescent 
eigenmodes using the conventional (o{k) method for calculating dispersion curves. 

4. Example: Negative Index Plasmonic Fishnet 

The second example highlighting the versatility of the CWES method is a plasmonic nega- 
tive index metamaterial (NIM) shown in Fig. [3] Because this metamaterial contains dispersive 
(plasmonic) components, it is even more convenient to use this method because the dielec- 
tric permittivities of metals are tabulated for real-valued frequencies. Recent experiments lfT9l 
demonstrated that the so-called fishnet metamaterial supports a negative index eigenmode for 
near infrared wavelengths of about Aq « 1.7 jxin. The dimensions and composition of the unit 
cell taken from Ref. [Wl are shown in Fig. |3] The fishnet metamaterial is made of alternating 
layers of Ag and MgFo with thicknesses of 3Qnm and 5Qnm respectively. This layered structure 
was milled with a focused ion beam into crisscrossing strips with widths of 265nm and 565«m. 
The crystal lattice constants are a, = = 860«m and a- = 80«m. For the dielectric response 
of the Ag we used the Drude model e^g = 1 — (0^/{(o{(0 — iy)) with the same parameters cited 
in Ref. = 9eV and y = 0.054ey. It should be noted that in Ref. the scattering 

frequency y was increased by a factor of three from that of standard bulk Ag to account for 
additional scattering losses. For the dielectric response of MgF2, no values were provided in 
the original paper so we used an oscillator model of the dielectric permittivity from Ref. ||251 . 

The negative index mode reported in Ref. |[T9| propagates in the z direction. Therefore, we 
have applied CWES to the specific case of ko = and k„ = z. The resulting dispersion curves 




Fig. 3. The fishnet metamaterial from Ref. 1191 . The lattice constants of the unit cell are 
Ux = ay = 860«ra and a, = SOnm. The fishnet is made up of alternating layers of Ag with a 
thickness of 30nm and MgF2 with a thickness of 50ivn. The widths of the the crisscrossing 
fishnet strips are bi = 265nm and b2 = 565nm. 

are plotted in Fig.|4]for the IQQTHz < (o/2n < 3Q0THz frequency range. For clarity we have 
only plotted the four eigenmodes that have the smallest values of lm{k^) and therefore the 
longest propagation lengths. The negative index mode found in Ref. |fT9l is labeled to indi- 
cate that it is transversely polarized with the electric field pointing in the x direction. We can 
see that it is a negative index mode because (according to the convention where fields are Bloch 
periodic with the exponential factor exp[i(a)f — k • x)]) a negative index mode should have a 
wavenumber whose real and imaginary parts have the same sign. In Fig. HI that describes the 
Ex mode for frequencies below lOOTHz- In addition, we observe three other negative index 
modes, though they generally have a shorter propagation lengths than the E^- mode. Of these 
three modes, the two labeled Ey are transversely polarized with the electric field pointing in the 
y direction and the mode labeled E- is longitudinally polarized with the electric field pointing 
in the z direction. 

The figure of merit (FOM) is a number commonly used to quantify the quaUty of a negative 
index material. It is often defined [61 as the ratio between the real and imaginary parts of the 
index of refraction. For eigenmodes of the fishnet crystal with real (O and a complex k point- 
ing in the z direction this is equivalent to the ratio between the real and imaginary parts of 
or FOM= Re(A:^)/Im(A:-). Without knowledge of the complex wavenumber of an eigenmode 
of a metamaterial, the figure of merit must be calculated indirectly. For example, in Ref. |fT9l 
the FOM is estimated numerically from simulated transmission through a fishnet sample. With 
knowledge of the complex wavenumbers of the crystal eigenmodes, it is possible to calculate 
the FOM directly. The FOM for each of the four eigenmodes identified in Fig. |4]is plotted in 
Fig-ISa), along with the propagation lengths plotted in Fig.|5jb). A positive FOM corresponds 
to the real and imaginary parts of k- having the same sign, therefore indicating a negative in- 
dex eigenmode. We see again that all four modes are in fact negative index modes in a fairly 
broad ISOTHz < (o/2n < IQQTHz frequency range. Somewhat unexpectedly, the mode exper- 
imentally observed in Ref. fT9l (labelled Ex) has the lowest FOM despite having the longest 
propagation length. The theoretically predicted FOM« 7.5 of the Ex mode is a factor of 2 
higher than the experimentally measured FOM |[T91 |, which could be attributed to fabrication 
imperfections. 

The existence of multiple negative index waves in the fishnet structure has important ex- 
perimental implications. For some frequencies (e.g., at ISQTHz) the propagation length of the 
longitudinal mode (£-) is as long as that of the primary transverse mode (Ex). While the orig- 
inal experiments ||T9l only excited the Ex mode by using the light normally incident onto the 




Fig. 4. TOP: Real (a) and imaginary (b) parts of k~{co) for several eigenmodes of the fishnet 
metamaterial shown in Fig.[3] In addition to the transverse mode electrically polarized in the 
X direction labeled which was identified in Ref. {19] we see two other transverse modes 
electricaUy polarized in the y direction labeled Ey and a longitudinal mode electrically 
polarized in the z direction labeled E,. MIDDLE: Field profiles for the transverse Ex mode 
(c) and the longitudinal E- mode (d) on a cross-section laying on the x-y plane in the middle 
of the MgF2 layer Arrows represent in-plane electric field and color represents the E, field. 
BOTTOM: (e) Field profile of the transverse E.r mode on a cross-section laying on the x-z 
plane halfway between two thin Ag-MgF2 strips. Arrows represent in-plane electric field 
and color represents the H,, field. For all field profiles the frequency is \15THz and each 
region is labeled Ag or MgF2 according the the material of the region. Unlabeled regions 
are vacuum. 



x — y vacuum/fishnet interface, one can envision exciting both and E- modes using obliquely 
incident light. For example, if p-polarized light is obliquely incident with the incident wavevec- 
tor laying in the x-z plane (electric field laying in the x-z plane), then both E^ and E^ modes 
will be launched into the fishnet with different phase velocities. Their refraction by the fish- 
net prism \\\9\ would give rise to two distinct beams. The existence of the additional strongly 
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Fig. 5. The POM (a) and propagation length (b) for the four fishnet eigenmodes identified 
in Fig. |4] Note that though the mode labeled has the largest propagation length it also 
has the smallest FOM. The negative value of FOM for the Ex mode for frequencies above 
llQTHz indicates that it is no longer a negative index mode at these higher frequencies. 

dispersive longitudinal modes (bulk plasmons) also suggests that the fishnet is a metamaterial 
with significant spatial dispersion that can have a strong effect on the surface waves at the vac- 
uum/fishnet interface 1261 [27l . Strong spatial dispersion is not unexpected because the fishnet's 
transverse period is about A /2. 

Recent theoretical ||28l and experimental ||29l studies demonstrated that negative index prop- 
agation in fishnet structures is not limited to the optical frequency range and can, in fact, be 
observed with microwaves. The main physical distinction (excluding their appropriately scaled 
sizes) between microwave and optical structures is that the metals in the former can be accu- 
rately described as PECs. On the other hand, in the optical range metals are described as being 
"plasmonic", which means that optical field penetration into the metal is significant. Our sim- 
ulation can quantify how important the plasmonic properties of the metal are for the fishnet 
structure from Ref. [T9l. This is done by introducing the plasmonic parameter Tp |f30 |, which 
is the number that characterizes the plasmonic nature of a metamaterial and it is defined as 
the ratio of the kinetic energy of the plasmonic electrons to the magnetic energy in the unit 
cell of a crystal. Strong plasmonic effects and the importance of electrostatic resonances imply 

^ 1. We find that the plasmonic parameter for the mode labelled Ej in Fig. 2] at a fre- 
quency of llSTHz is Tp = 0.24. Being less than one indicates this mode is not predominately 
plasmonic in nature because the electromagnetic fields do not significantly penetrate into the 
Ag. This is consistent with a recent demonstration of a negative index mode in a PEC fishnet 
structure ||29ll28l . 

Finally, we describe another application of the CWES method: calculation of isofrequency 
e)(k) = const contours in metamaterial crystals. From the early days of metamaterials reseach, 
isofrequency diagrams have been used as a simple visual tool for studying refraction at the vac- 
uum/metamaterial interfaces, especially in the context of negative index propagation and nega- 
tive refraction 1311 . Traditionally, isofrequency diagrams are drawn using a conventional ©(k) 
eigenvalue simulation fS). This provides no information about propagation lengths which are 
important in lossy metamaterials such as the fishnet. The conventional approach is also highly 
laborious for plasmonic fishnets because of the strong frequency dependence of the dielectric 
permittivities of metals. Other semi-analytic techniques used for analyzing wave propagation 
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Fig. 6. (a) Isofrequency contours for vacuum plotted with respect to k~ and ky (b) Isofre- 
quency contours for the E.v polarized eigenmode of the fishnet crystal plotted with respect 
to Re(/:j) and ky. (c) Isofrequency contours for the E j polarized eigenmode of the fishnet 
crystal plotted with respect to Im(/:.) and fc,,. The black arrows indicate the direction of the 
phase velocity and the red arrows indicate the direction of the group velocity at a frequency 
of nOTHz and ky = (o/cs'm{n/6). 

in highly symmetric PEC-based fishnets ll28l do not apply to plasmonic fishnets of interest. 

Figure |6] shows an isofrequency diagram calculated using the CWES approach. A single 
isofrequency contour was obtain by fixing the real frequency CO, setting ko = kyf and k„ = z, 
and then scanning ky from —n/oy to n/oy. The eigenvalue in this simulation is a complex 
valued k-. This procedure was repeated for several values of co from the l5QTHz to IWTHz- 
The resulting eigenmodes can be be excited by a plane wave incident on an interface between 
vacuum and the fishnet crystal parallel with the x-y plane if the wavevector of the incident 
wave is confined to the y-z plane. Because the y component of the incident wavevector is real 
valued, the y component of the wave excited in the fishnet crystal must also be real-valued. 

= sin 9 CO /c where 9 is the incidence angle with respect to the normal z. The eigenmode 
excited in the fishnet crystal decays in the z direction as indicated by the imaginary part of k^ 
plotted in Fig.|6] 

The isofrequency contours in Fig.|6]are hyperbolic in appearance. We can study refraction at 
the interface between vacuum and the fishnet crystal by comparing the isofrequency contours 
of the fishnet to those of vacuum, which are also shown in Fig. |6] The vacuum isofrequency 
contours are circular. As can be seen in Fig. |6]for frequency of llOTHz and with an incident 
angle of 30° the component of the incident wavevector tangential to the interface (ky) must be 
matched to the eigenmodes of the fishnet. For ky = co/csin{n/6) and a frequency of llOTHz 
the isofrequency diagram shows two eigenmodes with the correct value of ky and equal and op- 
posite values of k-. We find the correct mode by calculating the group velocity Vg = dco/dRe{k) 
which is by definition normal to the isofrequency contours. In the absence of anomalous dis- 
persion the group velocity indicates the direction of energy flow in the crystal ll32l . Choosing 
the correct solution requires us to choose the solution that has energy flowing in the positive z 
direction (i.e. away from the interface). This selects the solution with a negative Re(A:-). This is 
a negative index mode in the sense that in the z direction the phase velocity and group velocity 
have opposite signs. 

However, because the shape of the isofrequency contours is hyperbolic, the phase and group 
velocities in the y direction have the same sign. Therefore, positive refraction at the x-y vac- 
uum/fishnet interface is expected according to Fig.|6]despite the fact that a negative index eigen- 
mode is excited. This does not contradict the experiment by Valentine et al. in Ref. Iil9j . where 



the fishnet structure was cut in the shape of a prism. In that experiment, the wave propagation 
through the fishnet was entirely in the z direction enabhng the group and phase velocities to be 
antiparallel. Negative refraction indeed occurs at a vacuum-fishnet interface tilted with respect 
to the principal axis of the crystal. 

5. Conclusion 

In conclusion, we have presented a three-dimensional realization of the Complex Wavenum- 
ber Eigenvalue Simulation (CWES) approach to calculating photonic band structure of meta- 
material/photonic crystals. A detailed implementation of CWES using FEM discretization is 
described. The CWES approach was applied to two periodic photonic structures: (a) photonic 
crystal comprised of dielectric spheres, and (b) plasmonic fishnet metamaterial supporting neg- 
ative refractive index waves. For case (a) we have used the results of CWES to identify both 
transverse and longitudinal modes and investigated their coupling that gives rise to avoided 
crossings and bandgap formation. For case (b) we have identified, for the first time, four 
negative-index modes of the fishnet structure and computed hyperbolic isofrequency surfaces 
for the least-damped transverse mode. 
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